import ij.*;
import ij.process.*;
import ij.gui.*;
import java.awt.*;
import ij.plugin.frame.*;


import java.awt.event.*;
import java.lang.*;
import java.util.*;
import javax.swing.*;

public class Check_uniformity3_ extends PlugInFrame implements ActionListener,FocusListener 
{

/*version 2:
the principle of the program remains identical but:
check the snr versus "the mean of the snr values of 
all rois" instead of the "snr of the complete image" */

/*version 3: 
improved reporting*/

	JButton go;
	ImagePlus imp,impcut;
	ImageProcessor ip,ipres,ipcut,ipresmean,ipresstddev,ipresmin,ipresmax,ipresmode,ipressnr,ipresdevpix,ipresdevroi;
	ImageStack istackres;
	JTextField TFfragsize,TFdevroilim,TFdevpixlim,TFoffset;
	int cutsize,devroilim=15,devpixlim=15;
	float offset=0.0f;
	JTextArea TAinfo;
	JButton copy;
	Checkbox checkdevpixmap,checkunimap;
	boolean booldevpixmap,boolunimap;
	
	public Check_uniformity3_() 
	{
		super("Check uniformity");
		cutsize=100;
		JPanel panel1=new JPanel();
		panel1.setLayout(new GridLayout(6,2));
		panel1.setBackground(SystemColor.control);
		
	//	TFfragsize.setToolTipText("field size of image fragment");
		panel1.add(new JLabel("ROI size in pixels:")); 
		TFfragsize = new JTextField(Integer.toString(cutsize));TFfragsize.addFocusListener(this);TFfragsize.addActionListener(this);panel1.add(TFfragsize);
		panel1.add(new JLabel("Limit for deviating ROI (mean and SNR %):")); 
		TFdevroilim = new JTextField(Integer.toString(devroilim));TFdevroilim.addFocusListener(this);TFdevroilim.addActionListener(this);panel1.add(TFdevroilim);
		panel1.add(new JLabel("Limit for deviating pixels (%):")); 
		TFdevpixlim = new JTextField(Integer.toString(devpixlim));TFdevpixlim.addFocusListener(this);TFdevpixlim.addActionListener(this);panel1.add(TFdevpixlim);
		panel1.add(new JLabel("Offset value:"));
		TFoffset=new JTextField(Float.toString(offset));TFoffset.addFocusListener(this);TFdevpixlim.addActionListener(this);panel1.add(TFoffset);
		
		checkunimap=new Checkbox("Generate uniformity maps");
		panel1.add(checkunimap);
		
		checkdevpixmap=new Checkbox("Generate a map with deviating pixels");
		panel1.add(checkdevpixmap);
		
		go = new JButton("Check uniformity");
		go.addActionListener(this);
		panel1.add(go);
		copy=new JButton("Copy all data to system clipboard");
		copy.addActionListener(this);
		panel1.add(copy);
		add(panel1,BorderLayout.NORTH);
		TAinfo=new JTextArea(20,20);
		TAinfo.setEditable(false);
		JScrollPane JSP=new JScrollPane(TAinfo);
		add(JSP,BorderLayout.SOUTH);
		
		
		pack();
		GUI.center(this);
		show();
	}
	
	void showAbout()
	{ 
		IJ.showMessage("About Check Uniformity...", "This plugin ..." );
	}

	public void copyData()
	{
			JTextArea TA2=new JTextArea(TAinfo.getText());
			TA2.selectAll();
			TA2.copy();
			return;
	}

	double parseDouble(String s) throws NumberFormatException 
	{
		Double d = new Double(s);
		return(d.doubleValue());
	}	
	
	public void focusLost(FocusEvent e) 
	{
		if (e.getSource() == TFfragsize) 
		{
			try 
			{	cutsize = (int)Integer.parseInt(TFfragsize.getText());}
    			catch (NumberFormatException nfe) 
    			{	cutsize=100;}
    			
    			if(!(cutsize>0.0))
    				cutsize=100;
    			
    			TFfragsize.setText(Integer.toString(cutsize));
    			return;
		}
		
		if (e.getSource() == TFdevroilim) 
		{
			try 
			{	devroilim = (int)Integer.parseInt(TFdevroilim.getText());}
    			catch (NumberFormatException nfe) 
    			{	devroilim=5;}
    			
    			if((!(devroilim>0))||devroilim>100)
    				devroilim=5;
    			TFdevroilim.setText(Integer.toString(devroilim));
    			return;
		}
		
		if (e.getSource() == TFdevpixlim) 
		{
			try 
			{	devpixlim = (int)Integer.parseInt(TFdevpixlim.getText());}
    			catch (NumberFormatException nfe) 
    			{	devpixlim=15;}
    			
    			if((!(devpixlim>0.0))||devpixlim>100)
    				devpixlim=15;
    			TFdevpixlim.setText(Integer.toString(devpixlim));
    			return;
		}
		if (e.getSource() == TFoffset) 
		{
			try 
			{	offset = (float)Double.parseDouble(TFoffset.getText());}
    			catch (NumberFormatException nfe) 
    			{	
	    			offset=0.0f;
    			}
    			
    			TFoffset.setText(Double.toString(offset));
    			return;
		}
		
	}

	public void actionPerformed(ActionEvent e) 
	{
		ImageStatistics imstat;
		imp = WindowManager.getCurrentImage();
		int nx,ny,w,h,i,j,x0,y0,nrroi=0,ndevroi=0,ndevpix=0;
		char devpix,devroi;
		byte [] devpixmap;
		
		float imamean,imasnr,imastd;
		
		
		if (e.getSource()==copy)
		{
			copyData();		
			return;
		
		}
	
		
	
		if (imp == null) 
		{
			IJ.error("There is no active image, Please open an image.");
			return;
		}
		
		
		if (e.getSource() == go)
		{
			
			booldevpixmap=checkdevpixmap.getState();
			
			boolunimap=checkunimap.getState();
				
			TAinfo.setText("");
			
			float drl=(float)devroilim/100.0f;
			float dpl=(float)devpixlim/100.0f;
			cutsize = (int)Integer.parseInt(TFfragsize.getText());
			
			TAinfo.append("Uniformity analysis of "+imp.getTitle()+" based on "+cutsize+"x"+cutsize+" ROIs"+"\n");
			TAinfo.append("\n");
			TAinfo.append("Selected ROI-size: "+cutsize+"x"+cutsize+" pixels\n");
			TAinfo.append("Selected pixel value offset: "+offset+"\n");
			TAinfo.append("\n");
			TAinfo.append("Statistics of the full image:\n");
			
			ip=imp.getProcessor();
			w=ip.getWidth();
			h=ip.getHeight();
			nx=(int)Math.floor((double)w/(cutsize/2)-1);
			ny=(int)Math.floor((double)h/(cutsize/2)-1);
			
			devpixmap=new byte[w*h];
			for(i=0;i<w*h;i++) devpixmap[i]=0;
			
			istackres=new ImageStack(nx,ny);

			float flmean[]=new float[nx*ny];
			float flstddev[]=new float[nx*ny];
			float flsnr[]=new float[nx*ny];
			float flmin[]=new float[nx*ny];
			float flmax[]=new float[nx*ny];
			float flmode[]=new float[nx*ny];
			float fldevroi[]=new float[nx*ny];
			float fldevpix[]=new float[nx*ny];
			
			int roi_x0[]=new int[nx*ny];
			int roi_y0[]=new int[nx*ny];
		//	int roi_x1[]=new int(nx*ny);
		//	int roi_y1[]=new int(nx*ny);
			
			
			ip.setRoi(0,0,ip.getWidth(),ip.getHeight());
			ipcut=ip.crop();
			impcut=new ImagePlus("ipcut",ipcut);
			imstat=impcut.getStatistics();
			imamean=(float)(imstat.mean);
			imastd=(float)(imstat.stdDev);
			imasnr=(float)((imamean-offset)/imastd);
			
			TAinfo.append("Average:\t"+(Math.round(imamean*1000.0)/1000.0)+"\n");
			TAinfo.append("StdDev:\t"+(Math.round(imastd*1000.0)/1000.0)+"\n");
			TAinfo.append("SNR:\t"+(Math.round(imasnr*1000.0)/1000.0)+"\n\n");
			TAinfo.append("Statistics of each individual ROI:\n");
			TAinfo.append("x0\ty0\tmean\tstddev\tSNR\tmax\tmin\tdevroi\tdevpix\n");
			for(i=0;i<nx;i++)
				for(j=0;j<ny;j++)
				{
					
					x0=(int)(  ((w-(double)cutsize)*i)/(nx-1.0) );
					y0=(int)(  ((h-(double)cutsize)*j)/(ny-1.0) );
					
					roi_x0[i+nx*j]=x0;
					roi_y0[i+nx*j]=y0;
					
					ip.setRoi(x0,y0,cutsize,cutsize);
					ipcut=ip.crop();
					impcut=new ImagePlus("ipcut",ipcut);
					
					imstat=impcut.getStatistics();
					flmean[i+nx*j]=(float)(imstat.mean);
					flstddev[i+nx*j]=(float)(imstat.stdDev);
					flsnr[i+nx*j]=(float)((imstat.mean-offset)/imstat.stdDev);
					flmin[i+nx*j]=(float)(imstat.min);
					flmax[i+nx*j]=(float)(imstat.max);
					flmode[i+nx*j]=(float)(imstat.dmode);
				// generate map of rois with deviating pixels (pixel value deviates more than ...% from mean in roi )
					
					float under_pix=offset+(flmean[i+nx*j]-offset)*(1.0f-dpl);
					float upper_pix=offset+(flmean[i+nx*j]-offset)*(1.0f+dpl);
							
				
					if(((flmax[i+nx*j])>upper_pix)
						||(flmin[i+nx*j]<under_pix))
					{	fldevpix[i+nx*j]=1.0f;
						{
							for(int iroi=0;iroi<cutsize;iroi++)
								for(int jroi=0;jroi<cutsize;jroi++)
								{
									float p=ip.getPixelValue(x0+iroi,y0+jroi);
									if((p<under_pix)||(p>upper_pix))
										if(devpixmap[x0+iroi+(y0+jroi)*w]<1)
										{
											devpixmap[x0+iroi+(y0+jroi)*w]=100;
											ndevpix++;
										}
								}
						}		
					
					}
					else 
					{	
						devpix='-';fldevpix[i+nx*j]=0.0f;
					}
					
					nrroi++;
				}
			
			
			ipresmean=new FloatProcessor(nx,ny,flmean,null);
			istackres.addSlice("mean",ipresmean);

			ipresstddev=new FloatProcessor(nx,ny,flstddev,null);
			istackres.addSlice("stdev",ipresstddev);

			ipressnr=new FloatProcessor(nx,ny,flsnr,null);
			istackres.addSlice("snr",ipressnr);

			impcut=new ImagePlus("ipcut",ipressnr);
			imstat=impcut.getStatistics();
			
			float snrallrois=(float)imstat.mean;
			
			float under_mean=offset+(1.0f-drl)*(imamean-offset);
			float upper_mean=offset+(1.0f+drl)*(imamean-offset);
			float under_snr=(1.0f-drl)*snrallrois;
			float upper_snr=(1.0f+drl)*snrallrois;
					
			
			for(i=0;i<nx;i++)
				for(j=0;j<ny;j++)
				{
					
					x0=roi_x0[i+nx*j];
					y0=roi_y0[i+nx*j];
					
					
				// generate map of rois with deviating pixels (pixel value deviates more than ...% from mean in roi )
					if(fldevpix[i+nx*j]==1.0f)
					{	
						devpix='*';
					}
					else 
					{	
						devpix='-';
					}
				//generate map of deviating rois (mean and snr of roi deviates more than ...% from mean over image)		
					
					if ((flmean[i+nx*j]<under_mean)||(flmean[i+nx*j]>upper_mean)||
						(flsnr[i+nx*j]<under_snr)||(flsnr[i+nx*j]>upper_snr))
					{	
						devroi='*';
						if(flmean[i+nx*j]<under_mean)
							fldevroi[i+nx*j]+=1.0f;
						if(flmean[i+nx*j]>upper_mean)
							fldevroi[i+nx*j]+=2.0f;
						if(flsnr[i+nx*j]<under_snr)
							fldevroi[i+nx*j]+=4.0f;
						if(flsnr[i+nx*j]>upper_snr)
							fldevroi[i+nx*j]+=8.0f;
						
						ndevroi++;
					}
					else 
					{
						devroi='-';fldevroi[i+nx*j]=0.0f;
					}
					
					
					TAinfo.append(x0+"\t"+y0+"\t"+
						(Math.round(flmean[i+nx*j]*1000.0)/1000.0)+"\t"+
						(Math.round(flstddev[i+nx*j]*1000.0)/1000.0)+"\t"+
						(Math.round(flsnr[i+nx*j]*1000.0)/1000.0)+"\t"+
						(Math.round(flmax[i+nx*j]))+"\t"+
						(Math.round(flmin[i+nx*j]))+"\t"+devroi+"\t"+devpix+"\n"
					);
					nrroi++;
			
				}
			TAinfo.append("\n");
			TAinfo.append("Mean SNR in "+nrroi+" ROI's: "+(Math.round(imstat.mean*1000.0)/1000.0)+" with stdev: "+(Math.round(imstat.stdDev*1000.0)/1000.0)+"\n");
			TAinfo.append("Number of deviating ROI's ("+(devroilim)+"%):\t"+ndevroi+"\n");
			TAinfo.append("Number of ROI's with deviating pixel(s) ("+(devpixlim)+"%):\t"+ndevpix+"\n\n");
			
			ipresmin=new FloatProcessor(nx,ny,flmin,null);
			istackres.addSlice("min",ipresmin);

			ipresmax=new FloatProcessor(nx,ny,flmax,null);
			istackres.addSlice("max",ipresmax);

			ipresmode=new FloatProcessor(nx,ny,flmode,null);
			istackres.addSlice("dmode",ipresmode);

			ipresdevroi=new FloatProcessor(nx,ny,fldevroi,null);
			istackres.addSlice("deviating ROI's",ipresdevroi);
			
			ipresdevpix=new FloatProcessor(nx,ny,fldevpix,null);
			istackres.addSlice("deviating pixels",ipresdevpix);
			
			if(booldevpixmap)
			{	ByteProcessor ipdevpixmap=new ByteProcessor(w,h,devpixmap,null);
				new ImagePlus ("deviating pixel map",ipdevpixmap).show();
			}
			
			if(ndevroi>0)
			{
				TAinfo.append("Deviating ROI's\naccepted ranges:\nmean: ["
				+(Math.round(under_mean*1000.0)/1000.0)+","+(Math.round(upper_mean*1000.0)/1000.0)+
				"] \nsnr: ["+(Math.round(under_snr*1000.0)/1000.0)+","+(Math.round(upper_snr*1000.0)/1000.0)+"]\n");
				TAinfo.append("x\ty\tmean\tSNR\n");
				for(i=0;i<nx;i++)
					for(j=0;j<ny;j++)
					{
						if(fldevroi[i+nx*j]>0.0f)
							TAinfo.append(""+roi_x0[i+nx*j]+"\t"+roi_y0[i+nx*j]+"\t"+
							(Math.round(flmean[i+nx*j]*1000.0)/1000.0)+"\t"+
							(Math.round(flsnr[i+nx*j]*1000.0)/1000.0)+"\n");
					}
			
					
			}
			
			if(ndevpix>0)
			{
				ndevpix=0;
				TAinfo.append("Deviating pixels\n");
				TAinfo.append("x\ty\tpixel value\n");
				for(j=0;j<h;j++)
					for(i=0;i<w;i++)
						if (devpixmap[i+w*j]>0)
						{	TAinfo.append(""+i+"\t"+j+"\t"+ip.getPixelValue(i,j)+"\n");
							ndevpix++;
						}
			}
			
			
			
			if(boolunimap)
				new ImagePlus("Uniformity maps of "+imp.getTitle(),istackres).show();
		}	 
	}

}
